% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replication "Deconstructing the Yield Curve"
% Crump and Gospodinov (2024)
% Date: 26-JUL-2024
% Code for application on expectations hypothesis (Figures 4 and 5)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; clc; close all;
addpath('functions');
rng(123456);

o.dataSet               = 'gsw';
o.freq                  = 'q';
o.startDate             = 197201;
o.endDate               = 202401;
p.N                     = 40;
p.CG_B_BC               = 1000;

B = 1000;
M = 'ROT';

disp(' '); disp('Loading Data...'); disp(' ');
% Load data sets
dataAll = loadBondData(p.N);
dataTag = ['dataAll.', o.dataSet, '.', o.freq];
for i = {'prc', 'yld', 'ret', 'xret', 'fwd', 'dret'}, eval([i{1}, 'Raw = ', dataTag, '.', i{1}, ';']);  end
data.mats = (1:p.N);
datesBonds = eval([dataTag, '.dates']);
data.period = eval([dataTag, '.period']);

datesAll = 100*kron((1901:2026)',ones(4,1))+kron(ones(126,1),(1:4)'); %!%
tmpIdxStart = find(datesAll==o.startDate);
tmpIdxEnd = find(datesAll==o.endDate);
datesAll = datesAll(tmpIdxStart:tmpIdxEnd);

idxBonds = ismember(datesBonds,datesAll);

data.prc = prcRaw(idxBonds,1:p.N);
data.yld = yldRaw(idxBonds,1:p.N);
data.fwd = fwdRaw(idxBonds,1:p.N);
data.dret = dretRaw(idxBonds,1:p.N);
data.ret = retRaw(idxBonds,1:p.N);
data.xret = xretRaw(idxBonds,1:p.N);
data.T = size(data.prc,1);

vec_nm = NaN(p.N*(p.N-1)/2,2);
bhatAll = NaN(p.N*(p.N-1)/2,1);

out = bhatEH(data.fwd,data.yld);
vec_nm = out.vec_nm;

bhatAllBs1 = NaN(B,p.N*(p.N-1)/2);
bhatAllBs2 = NaN(B,p.N);
bhatAllBs3 = NaN(B,p.N);
bhatAllBs4 = NaN(B,p.N*(p.N-1)/2);

mats = (1:p.N);

if strcmp(M,'ROT'),  M = ceil((data.T*p.N)^(2/5)); end

for b = 1:B
    if mod(b,50)==0, disp([num2str(b), ' of ', num2str(B), ' completed']); end
    
    % CG Bootstrap
    forw0 = data.fwd;
    Dret = forw0(1:end-1,2:end) - forw0(2:end,1:end-1);
    varFwd = var1BC(forw0(:,end),p.CG_B_BC,0);    
    Zfwd = varFwd.VBC(2:end,:);
    Zinit = [Zfwd, Dret];    
    bsFwd = NaN(data.T,p.N);
    Z = [Zinit; Zinit(1:M-1,:)];
    idxRand = randi(data.T-M,[ceil(data.T/M),1]);
    idxBootstrap = kron(idxRand,ones(M,1)) + kron(ones(ceil(data.T/M),1),(0:M-1)');
    idxBootstrap = idxBootstrap(1:data.T-1);
    bsZ = Z(idxBootstrap,:);    
    bsFwd(:,end) = filter(1 ,[1 -varFwd.PhiBC],...
           [0; varFwd.muBC+(bsZ(:, 1)-mean(bsZ(:, 1)))],forw0(1,end));
    bsDret = bsZ(:,2:p.N);
    %bsFwd = buildForwards(bsDret, forw0(1,:), bsFwd(:,end));
    bsFwd(1,:) = forw0(1,:);
    for n = p.N-1:-1:1
        for t = 2:data.T
            bsFwd(t,n) = bsFwd(t-1,n+1) - bsDret(t-1,n);        
        end
    end
    bsYld = cumsum(bsFwd,2)./cumsum(ones(size(bsFwd)),2); 
    bsPrc = -ones(size(bsYld,1),1) * mats .* bsYld;
    bsFwd = [bsYld(:,1), (bsPrc(:,1:end-1)-bsPrc(:,2:end))];
   
    out_CG = bhatEH(bsFwd,bsYld);
    bhatAllBs1(b,:) = out_CG.bhat1;
    bhatAllBs2(b,:) = out_CG.bhat2;
    bhatAllBs3(b,:) = out_CG.bhat3;
    bhatAllBs4(b,:,:) = out_CG.bhat4;

end

idx_m1per_opt = 4;
idx_m1per = find(vec_nm(:,2) == idx_m1per_opt);
if strcmp(o.freq,'q')
    idx_m1yr = find(vec_nm(:,2) == 4);
    idx_m5yr = find(vec_nm(:,2) == 20);
else    
    idx_m1yr = find(vec_nm(:,2) == 12);
    idx_m5yr = find(vec_nm(:,2) == 60);
end

figNameTag = ['_' num2str(o.startDate) '_' num2str(o.endDate) '_M' num2str(M) '_m1per' num2str(idx_m1per_opt)];
datesPlots = datetime(round(datesAll/100), 3*(datesAll - 100*round(datesAll/100)), 15);
save(['../output/EH', '_' num2str(o.startDate) '_' num2str(o.endDate),'.mat']);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Charts for Figure 5: Time series of EH regressions: LHS and RHS
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Test 1: CS1 (n,m)=(40,4)
n = 40;
m = 4;
figName = figure;
plot(datesPlots,[[data.yld(m+1:end,n-m) - data.yld(1:end-m,n); NaN(m,1)], (m/(n-m))*(data.yld(1:end,n) - data.yld(1:end,m))], 'LineWidth', 1.5);
set(gca,'fontsize', 16);
legend('LHS', 'RHS', 'Location', 'SouthEast');
print(figName, '-dpdf', ['../output/Test1_timeseries_n', num2str(n), '_m', num2str(m), figNameTag, '.pdf']);
% Test 1: CS1 (n,m)=(40,20)
m = 20;
plot(datesPlots,[[data.yld(m+1:end,n-m) - data.yld(1:end-m,n); NaN(m,1)], (m/(n-m))*(data.yld(1:end,n) - data.yld(1:end,m))], 'LineWidth', 1.5);
set(gca,'fontsize', 16);
print(figName, '-dpdf', ['../output/Test1_timeseries_n', num2str(n), '_m', num2str(m), figNameTag, '.pdf']);

% Test 2: FB
n = 20;
figName = figure;
plot(datesPlots, [[movmean(data.yld((n-1)+1:end,1),[n-1 0])- data.yld(1:end-(n-1),1); NaN(n-1,1)], data.yld(1:end,n) - data.yld(1:end,1)], 'LineWidth', 1.5);
set(gca,'fontsize', 16);
print(figName, '-dpdf', ['../output/Test2_timeseries_n', num2str(n), figNameTag, '.pdf']);
% Test 3
n = 20;
figName = figure;
plot(datesPlots,[[data.yld((n-1)+1:end,1)-data.yld(1:end-(n-1),1); NaN(n-1,1)], data.fwd(1:end,n) - data.yld(1:end,1)], 'LineWidth', 1.5);
set(gca,'fontsize', 16);
print(figName, '-dpdf', ['../output/Test3_timeseries_n', num2str(n), figNameTag, '.pdf']);
% Test 4
n = 40;
m = 4;
figName = figure;
plot(datesPlots,[[data.fwd(m+1:end,n-m)-data.fwd(1:end-m,1); NaN(m,1)], data.fwd(1:end,n)-data.fwd(1:end,1)], 'LineWidth', 1.5);
set(gca,'fontsize', 16);
print(figName, '-dpdf', ['../output/Test4_timeseries_n', num2str(n), '_m', num2str(m), figNameTag, '.pdf']);
%
m = 20;
figName = figure;
plot(datesPlots,[[data.fwd(m+1:end,n-m)-data.fwd(1:end-m,1); NaN(m,1)], data.fwd(1:end,n)-data.fwd(1:end,1)], 'LineWidth', 1.5);
set(gca,'fontsize', 16);
print(figName, '-dpdf', ['../output/Test4_timeseries_n', num2str(n), '_m', num2str(m), figNameTag, '.pdf']);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Charts for Figure 4: Bootstrapped OLS estimates
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get y-axis limits
plot(vec_nm(idx_m1per,1),bhatAllBs1(:,idx_m1per), 'Color', .7*[1,1,1]);
ylim_benchmark1 = get(gca, 'yLim');
plot(vec_nm(idx_m5yr,1),bhatAllBs1(:,idx_m5yr), 'Color', .7*[1,1,1]);
ylim_benchmark2 = get(gca, 'yLim');
close all;
ylim_benchmark = [min(ylim_benchmark1(1),ylim_benchmark2(1)) max(ylim_benchmark1(1),ylim_benchmark2(2))];

% Test 1: CS1 (m=4)
figName = figure;
plot(vec_nm(idx_m1per,1),bhatAllBs1(:,idx_m1per), 'Color', .7*[1,1,1]);
hold on
plot(vec_nm(idx_m1per,1),out.bhat1(idx_m1per), 'k', 'LineWidth', 1.5);
yline(1, 'k', 'LineWidth', 1.5, 'LineStyle', '--');
set(gca,'fontsize', 16);
ylim(ylim_benchmark);
xlabel('$n$', 'Interpreter', 'Latex')
print(figName, '-dpdf', ['../output/Test1_1per', figNameTag, '.pdf']);
% Test 1: CS1 (m=20)
figName = figure;
plot(vec_nm(idx_m5yr,1),bhatAllBs1(:,idx_m5yr), 'Color', .7*[1,1,1]);
hold on
plot(vec_nm(idx_m5yr,1),out.bhat1(idx_m5yr), 'k', 'LineWidth', 1.5);
yline(1, 'k', 'LineWidth', 1.5, 'LineStyle', '--');
set(gca,'fontsize', 16);
ylim(ylim_benchmark);
xlabel('$n$', 'Interpreter', 'Latex')
print(figName, '-dpdf', ['../output/Test1_5yr', figNameTag, '.pdf']);
% Test 2: FB
figName = figure;
plot(bhatAllBs2', 'Color', .7*[1,1,1]);
hold on
plot(out.bhat2, 'k', 'LineWidth', 1.5);
yline(1, 'k', 'LineWidth', 1.5, 'LineStyle', '--');
set(gca,'fontsize', 16);
ylim(ylim_benchmark);
xlabel('$n$', 'Interpreter', 'Latex')
print(figName, '-dpdf', ['../output/Test2', figNameTag, '.pdf']);
% Test 3: CS2
figName = figure;
plot(bhatAllBs3', 'Color', .7*[1,1,1]);
hold on
plot(out.bhat3, 'k', 'LineWidth', 1.5);
yline(1, 'k', 'LineWidth', 1.5, 'LineStyle', '--');
set(gca,'fontsize', 16);
ylim(ylim_benchmark);
xlabel('$n$', 'Interpreter', 'Latex')
print(figName, '-dpdf', ['../output/Test3', figNameTag, '.pdf']);
% Test 4: CG (m=4)
figName = figure;
plot(vec_nm(idx_m1per,1),bhatAllBs4(:,idx_m1per,1)', 'Color', .7*[1,1,1]);
hold on
plot(vec_nm(idx_m1per,1),out.bhat4(idx_m1per,1), 'k', 'LineWidth', 1.5);
yline(1, 'k', 'LineWidth', 1.5, 'LineStyle', '--');
set(gca,'fontsize', 16);
ylim(ylim_benchmark);
xlabel('$n$', 'Interpreter', 'Latex')
print(figName, '-dpdf', ['../output/Test4_1per', figNameTag, '.pdf']);
% Test 4: CG (m=20)
figName = figure;
plot(vec_nm(idx_m5yr,1), bhatAllBs4(:,idx_m5yr,1)', 'Color', .7*[1,1,1]);
hold on
plot(vec_nm(idx_m5yr,1), out.bhat4(idx_m5yr,1), 'k', 'LineWidth', 1.5);
yline(1, 'k', 'LineWidth', 1.5, 'LineStyle', '--');
set(gca,'fontsize', 16);
ylim(ylim_benchmark);
xlabel('$n$', 'Interpreter', 'Latex')
print(figName, '-dpdf', ['../output/Test4_5yr', figNameTag, '.pdf']);

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EH Function
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = bhatEH(fwd,yld)

    N = size(fwd,2);

    out.vec_nm = NaN(N*(N-1)/2,2);
    out.bhat1 = NaN(N*(N-1)/2,1);
    out.bhat2 = NaN(N,1);
    out.bhat3 = NaN(N,1);
    out.bhat4 = NaN(N*(N-1)/2,1);

    j = 1;
    for n = 1:N    
        if (n>1)
            % Test: CS1
            tmpy = movmean(yld((n-1)+1:end,1),[n-1 0])- yld(1:end-(n-1),1);
            tmpx = yld(1:end-(n-1),n) - yld(1:end-(n-1),1);
            tmp = regress(tmpy,[ones(size(tmpy)), tmpx]);            
            out.bhat2(n) = tmp(2);
            % Test: FB
            tmpy = yld((n-1)+1:end,1)- yld(1:end-(n-1),1);
            tmpx = fwd(1:end-(n-1),n) - yld(1:end-(n-1),1);
            tmp = regress(tmpy,[ones(size(tmpy)), tmpx]);            
            out.bhat3(n) = tmp(2);
        end
        for m = 1:N
            if (m<n)
                % Test: CS2
                tmpy = yld(m+1:end,n-m)- yld(1:end-m,n);
                tmpx = m*(yld(1:end-m,n) - yld(1:end-m,m))/(n-m);
                tmp = regress(tmpy,[ones(size(tmpy)), tmpx]);            
                out.bhat1(j) = tmp(2);
                % Test: CG
                tmpy = fwd(m+1:end,n-m)-fwd(1:end-m,1);
                tmpx = fwd(1:end-m,n)-fwd(1:end-m,1);
                tmp = regress(tmpy,[ones(size(tmpy)), tmpx]);            
                out.bhat4(j) = tmp(2);
                % Keep maturity index values
                out.vec_nm(j,:) = [n,m];
                j = j+1;
            end
        end
    end
    
end    

